/**
 * \copyright
 * Copyright (c) 2020, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

/* reaction package header file */

//#ifndef reaction_INC
//#define reactions_INC

#include <vector>
// using namespace std;

// CB2406 #ifdef OGS_FEM_CAP // CAP_REACT

/* Structure for exchange of reaction rates */
class REACT_CAP
{
private:
public:
    REACT_CAP(void);   // Initialize
    ~REACT_CAP(void);  // Clean up

    // Data
    int* rateflag;   /* flag used for determining if reaction are calculated */
    long nodenumber; /* number of nodes, on which reactions are calculated */
    bool flag_cap;   /* flag if *.cap file exists   DL 28,10,08*/
    std::string data_file;
    std::string data_format;
    int mass_num;
    std::string mass_type;
    std::string species_define;
    std::vector<int>
        pcs_mass_idx;  // index of mass_trans order-number in processes list,
                       // (some other process are not mass transport)
    std::vector<int> pcs_ospecies_idx;  // index of output species order-number
                                        // in processes list     21.01.2009
    std::vector<int>
        pcs_nspecies_idx;  // index of nlog-output species order-number in
                           // processes list 02.02.2009
    std::vector<int> nspecies_idx;  // index of nlog-output species order-number
                                    // in NLOG_AC list

    std::vector<int> phase_list;
    std::vector<char*> phase_post;

    std::vector<char*>
        species_list_name;  // From .cap list section, species name
    std::vector<int>
        species_list_phase;  // From .cap list section, phase flag for species
                             // 0-solid, 1-gas, 2-aqueous
    std::vector<long int> species_list_idx;  // From .cap list section, index of
                                             // species in ChemApp DATA file

    std::vector<char*> species_name;  // From Process, list of species name for
                                      // chemical reaction
    std::vector<int> species_phase;   // From Process, phase flag for species
                                      // 0-solid, 1-gas, 2-aqueous
    std::vector<long int>
        species_idx;  // From Process, index of species in ChemApp DATA file
    std::vector<int> species_mobil;  // From Process, process mobile
    std::vector<int>
        species_dormant;  // From KINETIC_SPECIES, 0-equilibrium 1-kinetic

    std::vector<char*> species_kin_name;    // DL 03.2010
    std::vector<int> species_kin_phase;     //
    std::vector<long int> species_kin_idx;  //

    std::vector<char*> species_relative_activity_name;  // DL 03.2013
    std::vector<int> species_relative_activity_idx;
    std::vector<long int> species_relative_activity_idx_pcs;
    std::vector<double> species_relative_activity, species_relative_activity_ia,
        species_relative_activity_calc;
    std::vector<int> species_relative_activity_state;

    std::vector<int> Node_All_SI;
    bool show_data_file;

    std::vector<std::string> pcs_rename0, pcs_rename0_pre;
    std::vector<std::vector<std::string> > pcs_rename1, pcs_rename1_pre;
    std::vector<std::vector<double> > pcs_rename_stoi, pcs_rename_stoi_pre;
    std::vector<int> pcs_rename_idx0, pcs_rename_idx0_pre;
    std::vector<std::vector<int> > pcs_rename_idx1, pcs_rename_idx1_pre;

    int current_node;               // for Kin Reaction calc
    double current_TT, current_PP;  // for Kin Reaction calc
    typedef struct
    {
        std::string type;
        double logK;
        double rate;
        std::vector<char*> species_name;
        std::vector<int> species_phase;
        std::vector<long int> species_idx;
        std::vector<double> species_stoi;
        // for Kinetic reaction
        bool is_Kin;
        int logK_type;
        double logK_kin;
        double theta;
        double eta;
        double SA;
        double PF;
        std::vector<double> EA;
        std::vector<double> K25;
        std::vector<std::string> ref_species;
        std::vector<int>
            ref_idx;  // position in mass trans list --> node_ac matrix
        std::vector<double> ref_expo;
    } reaction;
    reaction React;
    std::vector<reaction> Kin_Reactions, Kin_HKF_Reactions;
    std::vector<int>
        Kin_HKF_index;  // look up table between Reactions and HKF_Species

    typedef struct
    {
        // input
        std::string name0;  // name in Helgeson datafile
        // from datafile
        std::string name;
        std::string abbrev;  // abbreviation
        std::string scform;  // structural chemical formula
        std::string ecform;  // elemental chemical formula
        std::string ref[2];
        double charge;  // charge
        int type;       // phase flag
        // 0-minerals that do not undergo phase transitions
        // 1-minerals that undergo one phase transition
        // 2-minerals that undergo two phase transitions
        // 3-minerals that undergo three phase transitions
        // 4-gases
        // 5-aqueous species
        double param[9][4];  // parameters
        // calculating result
        double S;   // entropy
        double H;   // enthalpy
        double G;   // Gibbs free energy
        double V;   // volume
        double Cp;  // heat capacity
    } SpeciesData;

    std::vector<std::vector<double> > node_logK;  // node, reaction
    std::vector<std::vector<double> > node_ac;    // node, spec
    std::vector<SpeciesData> Kin_HKF_species;
    std::vector<std::vector<double> > node_HKF_logK;

    std::vector<char*> nlog_name;          // correspond name with pcs and nlog
    std::vector<char*> species_nlog_name;  // use .cap to define -log(ac) value
    std::vector<int> species_nlog_phase;
    std::vector<long int> species_nlog_idx;

    int pcs_redox;
    std::vector<char*> species_redox_name;
    std::vector<int> species_redox_phase;
    std::vector<long int> species_redox_idx;
    std::vector<long> redox_stoi;
    bool check_no_reaction_nodes; /* flag if CheckNoReactionNodes has been
                                     performed */
    //	double temperature; /* temperature of water at node as specified in
    // input file */ 	int heatflag;      /* if 1, heat transport is active and
    // passed to phreeqc */
    std::vector<std::string>
        ic_2_bc_geometry_type;  // Convert eq concentrations at this gemetry to
                                // BC there
    std::vector<std::string> ic_2_bc_geometry_name;
    std::vector<size_t> ic_2_bc_GeoID;

    std::vector<std::string> ic_2_bc_species;
    std::vector<double> ic_2_bc_species_value;
    /**
     * @param geo_obj object of class GEOObjects managing the geometric entities
     * @param unique_name unique name to access the geometric entities in
     * geo_obj
     */
    void ConvertIC2BC(const GEOLIB::GEOObjects& geo_obj,
                      const std::string& unique_name);

    std::ofstream warning_out;
    // Member functions

    /**
     * @param geo_obj object of class GEOObjects managing the geometric entities
     * @param unique_name unique name to access the geometric entities in
     * geo_obj
     */
    std::ios::pos_type Read(std::ifstream*, const GEOLIB::GEOObjects& geo_obj,
                            const std::string& unique_name);

    void CreateREACT(void);
    void InitChemApp(void);
    void ResetElement(void);
    void DefSpeciesInChemApp(void);
    void DefSpeciesInChemAppList(void);
    void DefSpeciesInChemAppAll(void);

    void UndefOutSpecies(void);

    void RecoverChemSystem(void);

    bool CmpSpeciesList(void);
    bool CmpSpeciesListDAT(void);
    bool CmpSpeciesListPCS(void);
    bool CmpSpeciesListNLOG(void);
    bool CmpSpeciesListREDOX(void);
    bool CmpSpeciesListKinReact(void);  // DL 03.2010
    bool CmpSpeciesListKIN(void);       // DL 03.2010

    bool CmpSpeciesListRelativeActivity(void);  // DL 04.2013

    bool SetKinHKFspecies(void);  // DL 11.2010
    bool LoadHKFparam(void);
    bool HKFcalc(double T, double P);

    void LoadMassPCS(void);

    void CompNamePhase(std::string, char*&, int&);

    int CheckNoReactionNodes(void);
    void LoopNodeReact(int, int);
    std::vector<double> KineticReact(std::vector<double>);  // DL 05.10.10
    void KinInit(void);
    void KinRate(void);
    std::vector<double> derivs(std::vector<double>);  // DL 05.10.10
    std::vector<double> ODE(std::vector<double>, std::vector<double>, int,
                            double);
    void KinParamUpdata(int, int, std::vector<double>& spvc);
    void KinParamUpdataHKF(int, double T, double P, int err);
    void LoopNodeReact_Liquid_Vapor(int, int);
    void LoopNodeReact_Liquid_Solid(int, int);

    void ExecuteReactionsChemApp(int, int);  // DL 28,10,08

    void ShowDataFile(void);
    void SetAllSolidAsDormant(void);
};
extern std::vector<REACT_CAP*> REACT_CAP_vec;
extern bool REACT_CAP_Read(std::string, const GEOLIB::GEOObjects& geo_obj,
                           const std::string& unique_name);

#ifdef OGS_FEM_CAP  // CAP_REACT
#endif
